clear;close all;filebrowser;workspace;warning('off','all');
h(1) = figure('units','normalized','WindowStyle','docked','outerposition',[0 0 1 1]);
h(4) = figure('units','normalized','WindowStyle','docked','outerposition',[0 0 1 1]);
h(2) = figure('units','normalized','WindowStyle','docked','outerposition',[0 0 1 1]);
h(3) = figure('units','normalized','WindowStyle','docked','outerposition',[0 0 1 1]);
hc = 1.23984193;    %[eV*um]
eV2cm1 = 8065.54429;    %[cm-1/eV]
VR = 1000;          % [V] repeller voltage
AVMI = 46.7;        %[mm/sqrt(meV/V)]
LabelFontSize = 14;
% KEKRb2 = [400  240 15 8]; %[cm-1]
% RVMI = [(10+20)/2 (10+13)/2 (2.8+5.6)/2 (1.2+2.8)/2];
% RVMI = [13 10.5 6.2 3.6];
% RVMI_err = [3 2.5 1 1];

%%%-----------Baysian fit for KRb2----------------
disp('==========KRb2=============');
% UVwave =       [285  296  300  300  330  335  340  338  335]; %[nm]
% RVMIKRb2 =     [3.68 3.2  3.95 3.87 2.05 0.87 0.29 1.06 1.14];
% RVMIKRb2_err = [0.37 0.31 0.23 0.61 0.66 0.2  0.11 0.13 0.06];
% % RVMI = AVMI.*sqrt(KEKRb2./eV2cm1*1000/VR);
% EUV = hc./UVwave*1e3*eV2cm1;    %[cm-1]
% xdata = EUV;
% ydata = RVMIKRb2.^2;
% ydata_err = 2.*RVMIKRb2.*RVMIKRb2_err;
% fun =@(b,x) b(2).*(x-b(1));
% % Define default initial guess
% b0 = [3e4, (max(ydata)-min(ydata))/(max(xdata)-min(xdata))];
% % ydata_Err = (ydata0_Err*1e-11).*ydata.^2;
% W = 1./ydata_err.^2;
% [b,residual,jacobian,~,~,ErrorModelInfo] = nlinfit(xdata,ydata,fun,b0,'Weights',W);
% plot(EUV, RVMI.^2, 'ok', 'MarkerSize',7,'MarkerEdgeColor',[0.5,0.5,0.5],'MarkerFaceColor','r')
figure(h(1));
% UVwave =       [285  296  300  300  330  335 ...
%                 340  346  350  352  348  354  344 ...
%                 338  356  335]; %[nm]
% RVMIKRb2 =     [3.68 3.2  3.95 3.87 2.05 0.87 ...
%                 0.29 0.19 0.36 0.17 0.44 0.22 0.54 ...
%                 1.06 0.16 1.14];
% RVMIKRb2_err = [0.37 0.31 0.23 0.61 0.66 0.2 ...
%                 0.11 0.07 0.2  0.03 0.04 0.02 0.1  ...
%                 0.13 0.04 0.06];

%%-----Using Bayesian Analysis v2 to fit R^2------
UVwave =        [285    296    300    330  ...   % 10 Hz data
                    335    338    344    346 ...
                    348    352    354    340   342]; %[nm]
R2VMIKRb2 =     [14.182 11.56  13.733 4.554   ...
                    1.315  1.406  0.291  0.038...
                    0.193  0.031  0.043  0.134 0.445];
R2VMIKRb2_err = [2.559  2.462  3.39   2.089  ...
                    0.132  0.288  0.091  0.027...
                    0.030  0.011  0.007  0.066 0.107];
            
EUV = hc./UVwave*1e3*eV2cm1;    %[cm-1]
errorbar(EUV, R2VMIKRb2, R2VMIKRb2_err, ...
    'ok', 'MarkerSize',7,'MarkerEdgeColor',[0.5,0.5,0.5], ...
    'MarkerFaceColor','m', 'color','k', 'CapSize', 1, 'LineWidth', 0.5);
xlabel('Ionization photon energy (cm^{-1})');
xlim([2.8e4, 3.1e4]);
ylim([-1, 21]);
xlim([2.8e4, 3.6e4]);
ylabel('R^2 (mm^2)');
xt = get(gca, 'XTick');
set(gca, 'FontSize', LabelFontSize, 'XMinorTick', 'on');
xt = get(gca, 'YTick');
set(gca, 'FontSize', LabelFontSize, 'YMinorTick', 'on');

%%%-----fit the threshold---------
UVwavecutoff = 342;
xdata = EUV(UVwave <= UVwavecutoff);
ydata = R2VMIKRb2(UVwave <= UVwavecutoff);
ydata_err = R2VMIKRb2_err(UVwave <= UVwavecutoff);
fun =@(b,x) b(2).*(x-b(1));
% Define default initial guess
b0 = [3e4, (max(ydata)-min(ydata))/(max(xdata)-min(xdata))];
% ydata_Err = (ydata0_Err*1e-11).*ydata.^2;
W = 1./ydata_err.^2;
[b,residual,jacobian,~,~,ErrorModelInfo] = nlinfit(xdata,ydata,fun,b0);
% Calculate 95% confidence interval of the fitted parameters
ci = nlparci(b,residual,'Jacobian',jacobian);
disp(['threshold energy = ', num2str(b(1)), ' cm^{-1} or ', num2str(hc*1e3*eV2cm1/b(1)) , ' nm']);
disp(['Upper 95% Confident level', ' + ', num2str(ci(1,2)-b(1)), ' cm^{-1} or ', num2str(hc*1e3*eV2cm1/ci(1,2)) , ' nm']);
disp(['Lower 95% Confident level', ' - ', num2str(b(1)-ci(1,1)), ' cm^{-1} or ', num2str(hc*1e3*eV2cm1/ci(1,1)) , ' nm']);
xcurve = 2.9e4:(max(xdata)-min(xdata))/100:3.6e4;
xcurve1 = 2.8e4:(max(xdata)-min(xdata))/100:3.6e4;
ycurve = fun(b,xcurve);
hold on
plot(xcurve, ycurve, '-m', 'LineWidth', 1.5);
plot(xcurve1, zeros(size(xcurve1)), '--k', 'LineWidth', 1.5);

% %%%------plot R^2 in log scale------------
% figure(h(4));
% errorbar(EUV, RVMIKRb2.^2, 2.*RVMIKRb2.*RVMIKRb2_err, ...
%     'ok', 'MarkerSize',7,'MarkerEdgeColor',[0.5,0.5,0.5], ...
%     'MarkerFaceColor','m', 'color','k', 'CapSize', 1, 'LineWidth', 0.5);
% hold on
% plot(xcurve, ycurve, '-m', 'LineWidth', 1.5);
% xlabel('Ionization photon energy (cm^{-1})');
% % xlim([2.e4, 4e4]);
% ylim([5e-3, 50]);
% % xlim([2.8e4, 3.6e4]);
% ylabel('R^2 (mm^2)');
% xt = get(gca, 'XTick');
% set(gca, 'FontSize', LabelFontSize, 'XMinorTick', 'on');
% xt = get(gca, 'YTick');
% set(gca, 'FontSize', LabelFontSize, 'YMinorTick', 'on');
% % set(gca,{'XScale','YScale'},{'log', 'log'});
% set(gca,'YScale','log');
% %%%-----------Baysian fit for K2Rb----------------
% % disp('===========K2Rb============');
% % UVwave = [285 296 300]; %[nm]
% % RVMIK2Rb = [5.57 4.29 4.72];
% % RVMIK2Rb_err = [0.75 0.8 0.87];
% % EUV = hc./UVwave*1e3*eV2cm1;    %[cm-1]
% % errorbar(EUV, RVMIK2Rb.^2, 2.*RVMIK2Rb.*RVMIK2Rb_err, 'ok', 'MarkerSize',7,'MarkerEdgeColor',[0.5,0.5,0.5],'MarkerFaceColor','b');
% % 
% % xdata = EUV;
% % ydata = RVMIK2Rb.^2;
% % ydata_err = 2.*RVMIK2Rb.*RVMIK2Rb_err;
% % fun =@(b,x) b(2).*(x-b(1));
% % % Define default initial guess
% % b0 = [3e4, (max(ydata)-min(ydata))/(max(xdata)-min(xdata))];
% % % ydata_Err = (ydata0_Err*1e-11).*ydata.^2;
% % W = 1./ydata_err.^2;
% % [b,residual,jacobian,~,~,ErrorModelInfo] = nlinfit(xdata,ydata,fun,b0,'Weights',W);
% % 
% % ci = nlparci(b,residual,'Jacobian',jacobian);
% % disp(['threshold energy = ', num2str(b(1)), ' cm^{-1}']);
% % disp(['Upper 95% Confident level', ' + ', num2str(ci(1,2)-b(1)), ' cm^{-1}']);
% % disp(['Lower 95% Confident level', ' - ', num2str(b(1)-ci(1,1)), ' cm^{-1}']);
% % 
% % ycurve = fun(b,xcurve);
% % hold on
% % plot(xcurve, ycurve, '-b', 'LineWidth', 1.5);
% % legend({'KRb_2^+', 'fit to KRb_2^+', '0', 'K_2Rb^+', 'fit to K_2Rb^+'},'Location','northwest');
% 
% %%%-----------Tetramer spectroscopy over UV wavelength----------------
% figure(h(2));
% UVwave_unsort     = [350  340  352  346  348  354  344  338  342 356  335   347]; %[nm]
% Ncycle_unsort     = [200  201  200  530  400  1011 250  1075 963 339  1273  1004];
% UVpower_unsort    = [245  250  165  215  240  114  224  182  246 59   99    220];   %[mW]
% UVwx_unsort       = [363  373  353  357  353  347  362  384  381 347  376   365];  %[um]
% UVwy_unsort       = [314  341  361  278  349  363  314  259  335 360  333   318];  %[um]
% UVpower_unsort    = UVpower_unsort./UVwx_unsort./UVwy_unsort*350^2;
% 
% NK_unsort         = [19.4 33.8 4.8  97.9 56.9 99.3 75.1 311  115 1.8  345   33.8]*1e3;
% NKbkg_unsort      = [26   42   30   83   79   162  52   216  319 53   244   188];
% 
% NRb_unsort        = [13.8 6.6  6.7  49.4 29.2 56.1 24.0 262  97  6    258   40.5]*1e3;
% NRbbkg_unsort     = [20   15   17   60   33   80   41   107  156 35   140   103];
% 
% NKRb_unsort       = [15.0 8.8  6.0  16.7 17.9 30.8 19.6 188  45  5.1  168   32.4]*1e3;
% NKRbbkg_unsort    = [22   9    16   41   31   91   32   69   143 27   117   79];
% 
% NK2_unsort        = [39   25   10   50   39   129  42   168  182 40   172   126];
% NK2bkg_unsort     = [12   23   13   57   35   103  31   114  167 28   122   94];
% 
% NRb2_unsort       = [228  26   19   84   57   324  47   197  241 55   272   180];
% NRb2bkg_unsort    = [14   11   14   33   31   80   30   81   115 22   66    64];
% 
% NK2Rb_unsort      = [12   16   15   46   37   72   26   90   118 17   152   73];
% NK2Rbbkg_unsort   = [16   11   11   38   20   86   22   78   120 23   96    67];
% 
% NKRb2_unsort      = [21   21   33   43   85   113  36   157  122 29   285   77];
% NKRb2bkg_unsort   = [15   13   9    23   46   62   21   52   91  18   81    62];
% 
% NK2Rb2_unsort     = [138  28   112  60   225  301  77   159  203 120  479   89];
% NK2Rb2bkg_unsort  = [15   20   5    35   29   66   36   98   170 18   110   97];
% 
% NKRb2norm_unsort  = [161  360  802  511  1180 932  105 258   789 243  840   ];
% NK2Rb2norm_unsort = [1274 143  3682 1193 2918 3980 247 212  3021 3225 1339];
% 
% [UVwave, UVwave_order] = sort(UVwave_unsort);
% Ncycle = Ncycle_unsort(UVwave_order);
% UVpower = UVpower_unsort(UVwave_order);
% 
% NK   =   NK_unsort(UVwave_order);
% NKbkg = NKbkg_unsort(UVwave_order);
% NRb  =   NRb_unsort(UVwave_order);
% NRbbkg = NRbbkg_unsort(UVwave_order);
% 
% NKRb =   NKRb_unsort(UVwave_order);
% NKRbbkg =   NKRbbkg_unsort(UVwave_order);
% NK2  =   NK2_unsort(UVwave_order);
% NK2bkg = NK2bkg_unsort(UVwave_order);
% NRb2 =   NRb2_unsort(UVwave_order);
% NRb2bkg = NRb2bkg_unsort(UVwave_order);
% NK2Rb = NK2Rb_unsort(UVwave_order);
% NK2Rbbkg = NK2Rbbkg_unsort(UVwave_order);
% NKRb2 = NKRb2_unsort(UVwave_order);
% NKRb2bkg = NKRb2bkg_unsort(UVwave_order);
% NK2Rb2 = NK2Rb2_unsort(UVwave_order);
% NK2Rb2bkg = NK2Rb2bkg_unsort(UVwave_order);
% % NKRb2norm = NKRb2norm_unsort(UVwave_order);
% % NK2Rb2norm = NK2Rb2norm_unsort(UVwave_order);
% 
% EUV = hc./UVwave*1e3*eV2cm1;    %[cm-1]
% markersize = 5;
% subplot(2,2,[1 2])
% plot(UVwave, (NK2-NK2bkg)./Ncycle*200./UVpower*245, '-or');
% hold on
% plot(UVwave, (NRb2-NRb2bkg)./Ncycle*200./UVpower*245, '-oc', 'MarkerSize', markersize);
% plot(UVwave, (NK2Rb-NK2Rbbkg)./Ncycle*200./UVpower*245, '-og', 'MarkerSize', markersize);
% plot(UVwave, (NKRb2-NKRb2bkg)./Ncycle*200./UVpower*245, '-om', 'MarkerSize', markersize);
% plot(UVwave, (NK2Rb2-NK2Rb2bkg)./Ncycle*200./UVpower*245, '-ob', 'MarkerSize', markersize, 'MarkerFaceColor', 'b');
% hold off
% legend('(N_{K_2}-N_{K_2bkg})',...
%     '(N_{Rb_2}-N_{Rb_2bkg})',...
%     '(N_{K_2Rb}-N_{K_2Rb bkg})',...
%     '(N_{KRb_2}-N_{KRb_2bkg})',...
%     '(N_{K_2Rb_2}-N_{K_2Rb_2bkg})',...
%     'Location', 'northeastoutside');
% xlabel('Ionization laser wavelength (nm)');
% ylabel('Ion (count)');
% xt = get(gca, 'XTick');
% set(gca, 'FontSize', LabelFontSize, 'XMinorTick', 'on');
% xt = get(gca, 'YTick');
% set(gca, 'FontSize', LabelFontSize, 'YMinorTick', 'on');
% xlim([334 357]);
% ylim([-20 270]);
% subplot(2,2,[3 4])
% plot(UVwave, (NK-NKbkg)./Ncycle*200./UVpower*245, '-ok', 'MarkerSize', markersize,'MarkerFaceColor', 'b');
% hold on;
% plot(UVwave, (NRb-NRbbkg)./Ncycle*200./UVpower*245, '-ok', 'MarkerSize', markersize,'MarkerFaceColor', 'r');
% plot(UVwave, (NKRb-NKRbbkg)./Ncycle*200./UVpower*245, '-ok', 'MarkerSize', markersize,'MarkerFaceColor', 'k');
% hold off;
% xlabel('UV wavelength (nm)');
% ylabel('Ion (counts)');
% legend('(N_{K}-N_{K bkg})',...
%     '(N_{Rb}-N_{Rb bkg})',...
%     '(N_{KRb}-N_{KRb bkg})',...
%     'Location', 'northeastoutside');
% xlim([334 357]);
% ylim([0 1.5e5]);
% xt = get(gca, 'XTick');
% set(gca, 'FontSize', LabelFontSize, 'XMinorTick', 'on');
% xt = get(gca, 'YTick');
% set(gca, 'FontSize', LabelFontSize, 'YMinorTick', 'on');
% disp('=======================');
% 
% figure(h(3));
% markersize = 8;
% plot(UVwave, (NK2Rb2-NK2Rb2bkg)./Ncycle*200./UVpower*245, '-ob', 'MarkerSize', markersize, 'MarkerFaceColor', 'b');
% hold on
% plot(UVwave, (NKRb2-NKRb2bkg)./Ncycle*200./UVpower*245, '-om', 'MarkerSize', markersize);
% plot(UVwave, (NK2Rb-NK2Rbbkg)./Ncycle*200./UVpower*245, '-og', 'MarkerSize', markersize);
% % xpl=345.8:0.0001:346.2;
% % ypl = (xpl-346).*1e6;
% % plot(xpl,ypl,'-k', 'LineWidth', 5);
% hold off
% % legend('K_2Rb_2^+', 'KRb_2^+', 'K_2Rb^+', 'Location', 'northwest');
% xlabel('Ionization laser wavelength (nm)');
% ylabel('Ion (count)');
% % title('Ionization Spectroscopy of K_2Rb_2^*');
% xlim([334 357]);
% ylim([-20 270]);
% xt = get(gca, 'XTick');
% set(gca, 'FontSize', LabelFontSize, 'XMinorTick', 'on');
% xt = get(gca, 'YTick');
% set(gca, 'FontSize', LabelFontSize, 'YMinorTick', 'on');
% figure(h(2));
% savefig(h, ['plots/summaryplot_K2Rb2_spectroscopy.fig']);
saveas(h(1), ['plots/summaryBaysian_plot.png']);
% saveas(h(2), ['plots/summaryplot_K2Rb2_spectroscopy_Full.png']);
% saveas(h(3), ['plots/summaryplot_K2Rb2_spectroscopy.png']);
